package de.unibi.comet.examples;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import de.unibi.comet.fa.Alphabet;
import de.unibi.comet.fa.CDFA;
import de.unibi.comet.fa.DFAFactory;
import de.unibi.comet.fa.GeneralizedString;
import de.unibi.comet.fa.MarkovAdditiveChain;
import de.unibi.comet.util.Log;

public class pvalue {

	/** Determines the empiric character distribution of a given string. */
	public static double[] getEmpiricDistribution(Alphabet alphabet, String s) {
		Log.getInstance().startTimer();
		double[] freq = new double[alphabet.size()];
		for (int i=0; i<s.length(); ++i) {
			freq[alphabet.getIndex(s.charAt(i))]+=1;
		}
		for (int i=0; i<freq.length; ++i) {
			freq[i]/=(double)s.length();
		}
		Log.getInstance().stopTimer("Determine empiric distribution");
		return freq;
	}
	
	public static void main(String[] args) {
		// first 100 nt of c. glutamicum genome
		String genome = "gtgagccagaactcatcttctttgctcgaaacctggcgccaagttgttgccgatctcacaactttgagccagcaagcggacagtggattcgacccattga"; 
		
		Log.getInstance().setTimingActive(true);
		Log.getInstance().setLogLevel(Log.Level.VERBOSE);

		// the alphabet of nucleotides
		Character[] chars = {'a','c','g','t'};
		Alphabet alphabet =  new Alphabet(Arrays.asList(chars));
		// calculate empiric distribution from letter frequencies
		double[] charDist = getEmpiricDistribution(alphabet, genome);
		
		// create generalized string for a motif and its complement ...
		GeneralizedString p1 =  new GeneralizedString(alphabet, "ga?c");
		GeneralizedString p2 =  new GeneralizedString(alphabet, "g?tc");
		// ... and put them into a list
		List<GeneralizedString> l = new ArrayList<GeneralizedString>(2);
		l.add(p1);
		l.add(p2);
		
		// build cdfa (=counting deterministic finite automaton) from motifs
		CDFA cdfa = DFAFactory.build(alphabet, l);
		// optional step: minimize number of states (speeds up calculations)
		cdfa = cdfa.minimize();

		// count number of matches
		int matches = cdfa.countMatchesDFA(genome);
		
		// create markov additive chain from cdfa
		MarkovAdditiveChain mac = cdfa.createMAC(charDist); 
		// setup initial distribution
		double[][] dist = new double[cdfa.getStateCount()][matches+1];
		dist[0][0]=1.0;
		mac.setDistribution(dist);
		
		// optional (slower by approx. factor 10):
		// convert to logarithmic domain, then all calculations are done logarithmically
		// (results in logarithmic distribution)
		// mac.convertToLogSpace();
		
		// main: compute evolution of markov chain
		// this steps takes O(#matches * |genome| * #dfa-states)
		mac.step(genome.length());
		// read off result
		double pvalue = mac.getMassProbability(matches);
		
		System.out.println(String.format("dfa states: %d, matches: %d, p-value: %e", cdfa.getStateCount(), matches, pvalue));
	}
}
